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ABSTRACT 

We study the turbulence induced in the dust layer of a protoplanetary disk based on the energetics 
of dust accretion due to gas drag. We estimate turbulence strength from the energy supplied by 
dust accretion, using the radial drift velocity of the dust particles in a laminar disk. Our estimate 
of the turbulence strength agrees with previous analytical and numerical research on the turbulence 
induced by Kelvin- Helmholtz and/or streaming instabilities for particles whose stopping time is less 
than the Keplerian time. For such small particles, the strongest turbulence is expected to occur when 

the dust-to-gas ratio of the disk is ~ (h g /r) ~ 10~ 2 , where C e ff s» 0.2 represents the energy 
supply efficiency to turbulence and h g /r ~ 5 x 1(U 2 is the aspect ratio of the gas disk. The maximum 
viscosity parameter is a max ~ C c gT s (h g /r) 2 ~ 10 T s , where T s (< 1) is the non-dimensional stopping 
time of the dust particles. Modification in the dust-to-gas ratio from the standard value, 10~ 2 , by 
any process, results in weaker turbulence and a thinner dust layer, and consequently may accelerate 
the growth process of the dust particles. 

Subject headings: accretion, accretion disks — planets and satellites: formation — protoplanetary 
disks 



1. INTRODUCTION 

The first step of planet formation in protoplanetary 
disks is the collisional growth of (sub-)micron-sized dust 
particles (or aggregates), and consequent sedimentation. 
Dust settling and the formation of a dust layer at the 
midplane of the disk play an important role in the sub- 
sequent planetesimal formation process. Enhancement 
of the particle density in the dust layer accelerates the 
collisional growth rate. If the density enhancement is 
high enough, planetesimals may form through the grav- 
itational instability of the dust layer (Goldreich & Ward 
1973; Sekiya 1983). However, the turbulent motion of 
the gas hinders the dust layer from thinning. Dust par- 
ticles are stirred by turbulence and diffuse to high alti- 
tudes from the midplane. Thus, turbulence is an imped- 
iment to planetesimal formation. Turbulence excited by 
magneto-rotational instability (MRI) is so strong that 
the dust layer cannot become thin enough to induce 
planetesimal formation through gravitational instability 
(Johansen & Klahr 2005; Fromang & Papaloizou 2006; 
Turner et al. 2006; Carballidoet al. 2006, 2011; Fromang 
& Nelson 2009). In addition, the turbulence- induced col- 
lisional velocity can be high enough to destroy dust ag- 
gregates (Carballido et al. 2008, 2010), though several 
mechanisms have been proposed to overcome this diffi- 
culty (Johansen et al. 2007, 2011; Lyra et al. 2008, 
2009). Even if the gas disk exhibits an initially laminar 
flow (such as is expected in the dead zone where the ion- 
ization fraction of the gas is too low to couple to the mag- 
netic field, as mentioned in Gammie 1996 and Sano et al. 
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2000), dust settling itself induces turbulence. The dust 
particles tend to rotate around the central star faster 
than the gas because the gas experiences a pressure gra- 
dient force acting in the opposite direction of the gravity 
of the star. When the dust layer thins and the dust-to- 
gas ratio in the layer approaches unity, the velocity dif- 
ference from that of the upper gas layer induces Kelvin- 
Helmholtz (KH) instability and excites turbulence in the 
dust layer (e.g., Goldreich & Ward 1973; Sekiya & Ishitsu 
2000, 2001; Garaud & Lin 2004). The velocity differ- 
ence between the dust particles and the gas inside the 
dust layer also induces streaming instability (Youdin & 
Goodman 2005; Youdin & Johansen 2007; Johansen & 
Youdin 2007; see Chiang & Youdin 2010 for review on 
the various instability in the dust layer). Laminar gas 
disks are considered excellent sites for planetesimal for- 
mation. Thus, it is important to clarify how turbulence 
induced in the dust layer diffuses the dust particles, ver- 
sus dust sedimentation. Several analytical and numerical 
studies have focused on this problem (Cuzzi et al. 1993; 
Champney et al. 1995; Sekiya 1998; Dobrovolskis et al. 
1999; Johansen et al. 2006; Michikoshi & Inutsuka 2006; 
Weidenschilling 2006, 2010; Bai & Stone 2010a, 2010b). 

Sekiya (1998, hereafter S98) analytically solved for the 
structure of the dust layer under turbulence induced by 
KH instability. To derive the density profile of the dust 
layer analytically, the author adopted several assump- 
tions. First, the dust particles were assumed to be small 
enough and were coupled so tightly to the gas that the 
dust and the gas could be treated as a single fluid. Sec- 
ond, the density structure was adjusted to keep the dust 
layer marginally unstable to KH instability. Third, the 
effects of the Coriolis force and of the Keplerian shear 
were neglected. These assumptions allowed the author to 
determine analytically the density structure of the dust 
layer and to discuss how far from the midplane the dust 
particles diffuse due to turbulence. The analysis in S98 
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provided a useful guide for subsequent numerical studies 
that investigated turbulence of the dust layer with more 
realistic assumptions. In reality, however, dust particles 
that are not tightly coupled to gas play an important 
role in inducing turbulence. Also, other kinds of insta- 
bility, such as streaming instability, may occur before KH 
instability sets in (Bai & Stone 2010a). Thus, in addi- 
tion to performing numerical simulations, it is desirable 
to include an analytic discussion on the structure of the 
dust layer using a more general set of assumptions than 
adopted in S98. 

In this paper, we revisit the analysis presented in S98 
from a different point of view. In S98, the density struc- 
ture of the dust layer was determined by the condition 
that the dust layer was marginally unstable to KH insta- 
bility; i.e., the Richardson number, Ri, had a constant 
critical value. In the stability condition on Ri, the free 
energy due to the shear velocity is compared with the en- 
ergy needed to lift the material against the vertical grav- 
ity (Chandrasekhar 1961). In the case of a dust layer, 
the free energy of the azimuthal velocity is compared 
with the gravitational potential in the vertical direction. 
We also discuss the energetics of the dust layer to deter- 
mine the density structure, but in this paper we focus on 
the gravitational potential in the radial direction. Sus- 
taining steady turbulence requires an energy supply. The 
source of the free energy for KH instability or streaming 
instability is the velocity difference between the dust and 
the gas. Because this free energy is consumed in the pro- 
cess of exciting turbulence, a velocity difference must be 
continually induced to provide steady turbulence. This 
velocity difference originates from the force balance in 
the radial direction between the stellar gravity, the cen- 
trifugal force, and the pressure gradient. In other words, 
the velocity difference is induced because the gas resides 
in a slightly shallower effective potential (including both 
the stellar gravity and the gas pressure), than the effec- 
tive potential for the dust. The dust loses its angular 
momentum due to gas drag and drifts towards the star, 
while the gas gains angular momentum and drifts out- 
ward. The total angular momentum remains constant. 
Because the effective potential for the dust is deeper, 
this process releases gravitational energy that can be a 
source of the free energy for turbulence. (The relation- 
ship between the angular momentum variation, AL, and 
the energy variation, AE, for circular orbiting material 
is AE — HAL, where the angular velocity, 51, is different 
for the gas and for the dust.) If the gas disk is initially in 
a state of laminar flow, i.e., if gas drag is the only process 
that exchanges angular momentum, then dust accretion 
is the source of the free energy used to induce turbulence. 
The accretion rate of the dust due to gas drag has been 
calculated in the literature (Nakagawa et al. 1986, here- 
after NSH86; Weidenschilling 2003; Youdin & Chiang 
2004). Using this dust accretion rate, and assuming that 
a certain fraction of the accretion energy is transferred to 
turbulence, it is straightforward to calculate the strength 
of the turbulence and to determine the density structure 
of the dust layer. In this paper, we show that an analy- 
sis of the energetics based on the gravitational potential 
in the radial direction results in a dust layer structure 
that is qualitatively equivalent to the dust layer struc- 
ture developed in the analysis presented in S98. Most 
of the important properties of the S98 model are repro- 



duced. Our analysis avoids some of the assumptions that 
were adopted in S98. Specifically, we do not assume tight 
coupling between the gas and the dust, nor a marginally 
unstable structure for KH instability. (The effects of the 
Coriolis force and of Keplerian shear, which are still ne- 
glected in this paper, are discussed by Ishitsu & Sekiya 
2002, 2003; Gomez & Ostriker 2005; Chiang 2008; Bar- 
ranco 2009; Lee et al. 2010.) Thus, the results of this 
paper may be applied to more general situations, includ- 
ing situations in which the dust is weakly coupled to the 
gas, or in which streaming instability acts as a source of 
turbulence. 

In <J5J we describe our model assumptions. In the 
energy release rate of the accreting dust is calculated. 
In <21 the strength of the turbulence is derived and its 
asymptotic forms in the limits of a small dust-to-gas ra- 
tio and a tight dust-gas coupling are discussed. In fjH the 
results of our model are compared with the S98 model 
and with previous numerical simulations of KH instabil- 
ity and streaming instability. In ^6. 11 we discuss how 
dust layer formation reduces the radial drift velocity of 
the dust and the relative velocity of the dust particles. 
In §6.3[ we check the applicability of our model to turbu- 
lence in the dust layer for the KH instability and stream- 
ing instability cases. 

2. DISK MODEL AND ASSUMPTIONS 
2.1. Brief Model Description 

We consider a dust layer that forms after dust parti- 
cles have settled to the midplane of a gas disk around 
a star. Before the dust settles, the gas disk is assumed 
to be subject to laminar flow. This means that we con- 
sider the dead zone where MRI is inactive because of the 
low degree of ionization. Turbulent gas motion may be 
present even in the dead zone because sound waves prop- 
agate from the active layers at high altitudes and disturb 
the gas in the dead zone (Fleming & Stone 2003; Suzuki 
et al. 2010; Okuzumi & Hirose 2011). We assume that 
the dead zone is wide enough that the active layer cannot 
induce strong turbulence at the midplane. 

Even in an initially laminar disk, dust sedimentation 
and the formation of the midplane dust layer cannot pro- 
ceed in a perfectly undisturbed fashion. When the mid- 
plane dust-to-gas ratio reaches a critical value, the veloc- 
ity differences between the dust layer and the upper gas 
layer, or between the individual dust particles and the 
surrounding gas, will start to induce hydrodynamical in- 
stabilities, such as KH instability or streaming instability 
(see Chiang & Youdin 2010 for review). Consequently, 
turbulent diffusion of the dust particles terminates fur- 
ther dust settling. It is expected that a steady dust layer 
forms, in which particle settling and turbulent diffusion 
balance each other. In such a state, vertical settling is 
no longer an energy source for turbulence, but the radial 
accretion of the dust still proceeds. The turbulence is 
maintained by the energy liberated from dust accretion 
to the star. 

We calculate the dust accretion rate and consequent 
energy release, assuming that a steady dust layer has 
formed. The strength of the turbulence is then esti- 
mated. The balance between the turbulent diffusion and 
the dust settling determines the dust layer thickness. 
Thus, the structure of the dust layer, the dust accretion 
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rate, and the turbulence strength must be determined 
self-consistently. 

For simplicity, we focus on a narrow axisymmetric ring 
region at a certain radius r from the central star. We 
consider the dust layer structure only in the vertical di- 
rection, neglecting any radial variation in the properties 
of the gas disk and the dust layer. All the dust particles 
are assumed to have a uniform size, i.e., the size distri- 
bution of the dust particles is neglected in this paper. 

2.2. Vertical Structure of the Dust Layer 

If the effect of gas drag on the dust were negligible, 
the dust particles would orbit at the Keplerian velocity 
£!k = sj (GM)/r 3 , where G is the gravitational constant 
and M is the mass of the central star. However, the or- 
bital velocity of the gas is slightly less than the Keplerian 
velocity because the direction of the gas pressure gradi- 
ent is usually outward from the star. This cancels part of 
the contribution to the orbital velocity from the gravity 
of the star. Given a negligible influence of the drag force 
from the dust, the angular velocity of the gas would be 



GM 



(1 - 2rf) « n K (l - rf). 



(1) 



Here, r\ represents the deviation of the gas orbital velocity 
from the Keplerian value. It is half of the ratio between 
the gas pressure gradient force and the gravity of the 
central star, 

1 dP 

71 2p g rn 2 K dr ' [ ' 

where p g and P are the gas density and pressure, respec- 
tively. 

We consider a thin dust layer in which most of the 
dust particles have settled on the midplane of the gas 
disk. The gas disk is isothermal in the vertical direction 



with a scale height of h g 



,/Ok, where c s is the gas 



sound speed. Because the thickness of the dust layer is 
much smaller than the scale height of the gas disk, the 
gas density p g in the dust layer is treated as a constant. 
Its value is given by the column density of the gas disk 



as 



Pg 



'2izh„ 



(3) 



The dust density profile in the vertical direction, Pd(z), 
is determined by the balance between the dust settling 
and turbulent diffusion (Youdin & Lithwick 2007, here- 
after YL07) 0. The dust particles have a stopping time, 

Tstop = Ts^k > m which the velocity difference from the 
background gas flow becomes 1/e times due to gas drag, 
where T s is the non-dimensional stopping time. In a 
turbulent gas disk with a turbulent diffusion coefficient 
D g , the density profile of the dust layer composed of 

4 Note that the diffusion coefficient given by YL07 is derived 
assuming that the dust particles can be treated as passive parti- 
cles. Problems may arise if this formula is used in a case where 
the dust layer is so dense that its dust-to-gas ratio is larger than 
unity. Although a comparison of our model with the simulation by 
Johansen et al. (2006) shows good agreement even for a dust-to- 
gas ratio as large as 40 (see Fig. 0, modeling the effect of the dust 
inertia will be a subject for future investigation using a refinement 
of the present model. 



single-sized particles is Gaussian, pd — pd.fi ex P( — z 2 ), 
where z — z/(y2hd) is the non-dimensional vertical 
coordinate normalized by the dust scale height hd ~ 
[D g /(n K T s )} 1/2 . Note that the turbulent diffusion coef- 
ficient in the "z-direction" , D gi may have a significantly 
different value from the usual turbulent viscosity coeffi- 
cient i/ acc in accretion disks, where v &cc comes from the 
"r(9-component" of the Reynolds stress (e.g., Lesur & 
Ogilvie 2010). Nevertheless, we follow the conventional 
"a-prescription" , and express the diffusion coefficient as 
D g = aCghg for simplicity. This prescription and equa- 
tion (24) of YL07 give 



(4) 




(5) 



The dust density profile is written as 

Pd = fmidPg exp(-5 2 ) , 

where the midplane dust-to-gas ratio / m id is increased 
by a factor h g /hd from the total dust-to-gas ratio or the 
"metallicity" of the disk Zdi S k = S^/Sg, 



mid 



7 h 9 
^diskl — 



^disk 




(6) 



We define (3{z) as the ratio of the total (dust and gas) 
density to the gas density, 



Pg +Pd 
Pg 



1 + /mid exp( 



-5 2 ) 



(7) 



2.3. Typical Values of the Model Parameters 

In our model, the turbulent parameter a and the struc- 
ture of the dust layer (i.e., the thickness, hd, and the mid- 
plane dust-to-gas ratio, /mid) are determined for a given 
parameter set (the stopping time of the dust particles T s , 
the disk metallicity i?disk, and the parameters for the gas 
disk). The non-dimensional variables of the result (a, 
hd/r, and / m id) depend on the gas disk parameters only 
through 



(8) 



The numerical 



where fj ~ 77 ~ (h g /r) 2 ~ 10~ 3 — 10" 
calculations in the subsequent sections mainly use a value 
fj = (O.O5/7) 2 = 9 x 10 -4 , where 7 = 5/3, in order 
to compare our model with the simulation presented by 
Johansen et al. (2006, hereafter JHK06). In fx2~2l fj = 
0.05 2 is used for comparison with the results by Bai & 
Stone (2010), and in JEB fj = 2.92 x 10~ 3 is used to 
calculate the radial drift velocity of the dust in the gas 
disk model by Hayashi (1981). The fj dependence of the 
results will be discussed in 



3. TURBULENCE ENERGY SUPPLY 

The energy supply for turbulence comes from libera- 
tion of the gravitational energy of the dust. As the dust 
falls towards the central star, the dust particles pene- 
trate more deeply into the potential well of the star. 
Although the gas drifts outward to conserve the total 
angular momentum, the difference in the effective grav- 
itational potential between the dust (—GM/r) and the 
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gas (— GM(1 — 2rj)/r), including the work done by the 
pressure gradient, causes a net energy liberation. Part 
of the liberated energy is converted directly into ther- 
mal energy, and part is used for supplying energy to the 
turbulence (see discussion in i)6.3p . 

Dust particle accretion occurs either individually or 
collectively. Individual dust particles suffer gas drag 
given a velocity difference between the particle and the 
gas. Usually the orbital velocity of the gas is slower than 
that of the dust particle, and the gas drag force deceler- 
ates the particle orbital motion. Thus, the particle loses 
its angular momentum and drifts inward. In addition to 
this individual drift, a collective drag works on the whole 
dust layer. The dust layer at the midplane orbits faster 
than the gas layer at higher altitudes because in the dust 
layer the increased inertia of the enriched dust particles 
weakens the effect of the gas pressure. The dust layer 
rotates with a velocity close to the Keplerian value. The 
upper gas layer, devoid of the inertia of the dust, orbits 
more slowly than the midplane dust layer. If the gas 
exhibits turbulent viscosity, then the slower orbiting gas 
layer exerts a drag force on the dust layer, and conse- 
quently the dust layer loses its angular momentum and 
accretes inward. We consider the energy release rate of 
the dust caused by individual and collective drag. 

3.1. Dust Accretion Caused by Individual Drag 

We calculate the gravitational energy released when 
dust particles accrete towards the star due to drag on 
individual particles. In this calculation, we assume that 
the gas disk is in a laminar flow state, and that dust 
particles drift inward steadily. Any random motion of 
particles caused by gas turbulence is neglected. We ex- 
pect that in turbulent gas disks, the average motion of 
the dust particles can be estimated from the motion in 
the laminar disk (Bai & Stone 2010). Because the in- 
dividual drag is more effective than the collective drag 
when the midplane dust-to-gas ratio is less than unity 
(see Fig. Q] below) , we take the / m id <C 1 limit in the 
following calculation. The calculation for a general / m id 
is described in Appendix [Bj 

In a laminar disk, the particle drift velocity vd,r is given 
by equation (2.11) of NSH86. For pd <C p g , the particle 
radial velocity is (see also eq. [23] of Takeuchi & Lin 
(2002); note the factor of 2 difference in the definition of 

v) 

2T„ 



I'd,-, 



T? + l 



(9) 



The gas drifts in the opposite direction with a velocity 
v g<r . From angular momentum conservation, the gas drift 
velocity is 

Vg, r = ~ — Vd,r ■ (10) 

The effective gravities (including the pressure gradient 
force) acting on the dust and on the gas are gd = 
—GM/r 2 = -rfi| and g g = -GM(1 - 2 V )/r 2 = -rQ 2 g , 
respectively. While the total angular momentum is con- 
served, the total energy is not conserved because of the 
difference in the effective gravities acting on the gas and 

^ 0, providing a source 



on the dust, p d gdVd,r + P Q 
of energy for turbulence. 



J 9,r 



Note that, in the above calculation, Vdr and Ug^r are 
the "terminal velocities". This means that the gravi- 
tational accelerations of the dust and gas are balanced 
by the drag forces, and liberated energy translates di- 
rectly into thermal energy. Thus, no energy would be 
contributed to turbulence. In reality, however, if the gas 
disk is turbulent, a steady terminal velocity is not ex- 
pected, and acc elera tion phases of the dust must occur, 
as discussed in ^6.31 In an acceleration phase, the work 
done by gravity first provides kinetic energy, which can 
be translated into turbulent energy. We estimate this 
energy input to turbulence. In the following calculation, 
the energy liberated from the accreting dust is estimated 
using the terminal velocity described above for simplic- 
ity. Only a certain fraction of the liberated energy goes 
into turbulence. Thus, the energy input to turbulence is 
a factor C g(< 1) times the following estimation. The 
factor C e ff will be determined in 8 35.21 to be ~ 0.2 by 
comparing our model with the numerical simulation of 
turbulence in the dust layer developed by JHK06. 

Given that some portion of the accretion energy of the 
dust is consumed by the outward motion of the gas, the 
liberated gravitational energy per unit surface area of the 
disk is, 



dE dl 



dt 



{pdgdVd,r+P g g g V a ,r)dz 



(11) 



where the factor of 1/2 comes from the fact that half 
of the work done by gravity is used for acceleration (and 
deceleration) of the azimuthal velocity of the dust (and of 
the gas) as their semi-major axes change. Using equation 
((9]) and ([T0| . the energy liberation rate reduces to 



ggdra 

dt 



= 2n 2 vlVL K T s Y,d. 



drag 



(12) 



where the effective "surface density" Ed, drag of the dust 
is ^ 

£<i,drag = T 2 + i ^ d ■ ( 13 ) 

In the above calculation, we assume / m id <C 1. The 
calculation of the energy liberation rate for general / m id 
is described in Appendix|B] In the numerical calculations 
in the subsequent sections, we use equation (|B5j) for the 
effective surface density 0, 



-^drag 



exp(— z 2 ) 
T 2 + P 2 



T 2 



2{T 2 + P 2 



dz 



(14) 

If the dust particles are tightly coupled to the gas (T s <C 
1) and dust sedimentation is weak (/ m id -C 1), then 
£,i,drag i s simply equal to the dust surface density E^. 

3.2. Dust Accretion Caused by Collective Drag 

Next, we consider dust accretion due to collective drag 
acting on the entire dust layer. Because the faster- 
orbiting dust particles drag the gas in the dust layer, the 

5 Equation d 1 41) does not coincide with equation <l X 3 1) in the limit 
°f /mid <S 1- This disc repancy comes from the fact that the calcu- 
lation of equation I I14II includes corre ctio n terms of the order of r/ 2 
in the drift velocity, while equation d 1 3 1) con siders onl y the terms 
of r\. The difference between equations d 1 3 1) and d 141) is at most 
factor 2 (for T a ^> 1), and is not significant. 
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orbital velocity of the gas is largest at the midplane and 
decreases with altitude. If the gas in the dust layer is tur- 
bulent, the variation in orbital velocity v g ^g with altitude 
z induces Reynolds stress Pg z . This causes a transfer of 
angular momentum from the dust layer to the upper gas 
layer, resulting in accretion of the dust layer. 

To calculate the energy liberation rate, we make a 
few key assumptions. First, we consider only the 9z- 
component of the Reynolds stress, Pg z , neglecting the 
other components P r e and P zr . Ignoring P r g means that 
turbulence in the dust layer does not transfer the angu- 
lar momentum efficiently in the radial direction. This 
is expected for turbulence induced by hydrodynamical 
instabilities such as convective instability (e.g., Stone & 
Balbus 1996; Lesur & Ogilvie 2010). The Reynolds stress 
P zr is also neglected for simplicity. Brauer et al. (2007) 
pointed out that P zr changes the velocity profiles of the 
dust and gas from those derived by NSH86 by a factor 
~ 3, and thus P zr cannot be neglected in a rigorous dis- 
cussion. Obtaining the exact velocity profiles including 
P zr requires numerical calculations. In this paper, we use 
the velocity profiles calculated analytically by NSH86 for 
simplicity. This induces an error of a factor ~ 3 in our 
estimation of turbulence strength. Thus, our estimate is 
limited to an order-of-magnitude argument. 

The second assumption is that the turbulent layer has 
a thickness comparable to that of the dust layer. If the 
turbulent layer were much thicker than the dust layer 
and most of the volume of the turbulent layer were free 
of the dust, then its structure would be controlled by 
the gas, unaffected by the properties of the dust and the 
structure of the dust layer. The dust layer would behave 
just like a boundary wall at the bottom of the turbulent 
layer. Turbulence in such a thick boundary gas layer 
has been discussed using the analogy of the Ekman layer 
(e.g., Cuzzi et al. 1993). If the turbulent layer were 
dominated by the dust layer, then the structure of the 
dust layer would control the turbulence strength. We 
focus on such a dusty turbulent layer. Following Youdin 
& Chiang (2004), we consider the conditions needed for 
the turbulent and dust layers to be of similar thickness. 
The thickness of the turbulent layer, from dimensional 
analysis, is the Ekman length, 



(15) 



where v is the turbulent viscosity, provided that the vis- 
cosity and the Coriolis force determine the layer struc- 
ture. The thickness of the dust layer (eq. [4]), which is 
determined by the balance between sedimentation and 
diffusion of the dust particles, is 



T s Ok 



(16) 



where v ~ D g is used. For small particles (T s < 1), 
hd is larger than He, meaning that such small particles 
move further out of the turbulent layer and modify the 
structure of the turbulent layer. It is expected that the 
thickness of the turbulent layer is not given by He, but 
by hd (see also S98; Goodman & Pindor 2000). For large 
particles (T s ^> 1), the turbulent layer is much thicker 
than the dust layer, and its thickness is expected to be 
He- In the following discussion, we consider a dusty tur- 



bulent layer with a thickness similar to hd- Thus, the 
analysis in this paper is probably not appropriate for 
large particles (T s >• 1). 

Under the above assumptions, the energy liberation 
rate is estimated. Because the collective drag is effective 
only if the midplane dust-to-gas ratio is larger than unity 
(see Fig. rjjbclow), we consider the case in which a dense 
dust layer has formed (/ m id > 1) and use the plate drag 
approximation (Goldreich k. Ward 1973; Goodman & 
Pindor 2000; Weidenschilling 2003). The calculation for 
general / m id is described in Appendix [Cj The Reynolds 
stress Pg z near the boundary between the dust layer (or 
the turbulent layer) and the gas layer is estimated as 



Pez = P g v- 



dv 



y.O 



dz 



-PgV- 



hd 



(17) 



This stress extracts angular momentum from the dust 
layer and transfers it to the gas layer. The unit sur- 
face of the dust layer loses angular momentum dLd/dt = 
rPg z , and the corresponding energy change is dEd/dt = 
VlyidLd/ dt. The gas layer gains the same amount of 
angular momentum dL g /dt = —rPg Zl but the energy 
change is different from that of the dust layer because 
of the work done by the pressure gradient: dE g /dt — 
Vl g dL g /dt = (1 — r/jn^dLg/dt. In total, the energy lib- 
eration rate (the minus sign is added), using equations 
©,©,©, and (H7J), is 



dE vis 
dt 
where 



= -(fi K - Q g )rP ez = 2r] 2 v K 2 tt K T s Y, d s 



1 1 + 2T, 1 



2V2T 1 + T s f, 



-'d ■ 



mid 



(18) 



(19) 



In the above calculation, the viscosity coefficient is mod- 
eled as v — ac g h g . The same a is used for both D g 
and v for simplicity. The effective surface density £d, v is 
is inversely proportional to the midplane dust-to-gas ra- 
tio /mid and only weakly depends on T s . Note that the 
above relationship is derived for / m id > 1. The effective 
surface density for general / m id is derived in Appendix [C] 
and is given by 



^cLvis / — Jmid^-'d , m 
\/7T 1 + T s 



2 + 2/3T 2 - T 2 



z 2 exp(— 2z 2 ) 



(T s 2 + /3 2 ) 2 [1 + C str / mid exp(-z 2 )] 



where 



T? + l 



(20) 



(21) 



Expression (|20| is complicated, but for large / m id, its de- 
pendence on T s and / m id is similar to that of the simpler 
equation (fl~9| ; Ed, vis oc an< ^ depends only weakly 

on T s . 

3.3. Energy Dissipation in Turbulence 

The turbulent energy of the largest eddies transfers to 
smaller eddies, and finally dissipates to thermal energy 
via decay of the smallest eddies due to molecular viscos- 
ity. The size and velocity of the largest eddies are as- 
sumed to be ^ ff ,oddy = a 1 / 2 /i g and u g , c ddy = ot x l 2 c s . This 
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assumption means that the turnover time of the largest 
eddies is the Keplerian time (r Sj eddy = ^cddy/wg.cddy = 
; Cuzzi et al. 2001). Using the Kolmogorov scaling 
law, the energy dissipation rate per unit volume and unit 
time is 



di 



turb 



at 



u 



ff.cddy 
r g,cddy 



(22) 

where the factor C cnc represents the fact that dust par- 
ticles that are coupled only weakly to the gas do not 
contribute to the turbulent energy, and is given by (see 
Appendix [D]) 



C' 



T,T F 



T s 2 (Te 



for T s 
for T, 



(23) 



where T e = Tg ie ddy^K is the non-dimensional turnover 
time of the largest eddies. For most of this paper, we 
consider T e = 1. The effect of the dust inertia is ignored 
for simplicity. 

Energy dissipation of turbulence occurs only in the 
dust layer. This assumption may be problematic for 
T s S> 1, as discussed in ^3 . 2 1 but for simplicity, we as- 
sume that energy dissipation occurs in ~\/2hd < z < 
V2hd- The dissipation rate of the energy per unit sur- 
face area and unit time is 



<9£-turb ,2 n 3 



\/2h t 



V2h d 



(p g + C C ncPd)dz = ah 2 Vl\Y>d,tuA ) 



where 

Sd,turb = 

and erf (1) = 0.8427. 



V^/mid 



+ C ono erf (1) 



(24) 



(25) 



4. TURBULENCE STRENGTH 

The strength of turbulence or the parameter a is de- 
termined by the balance between the energy supply rate 
(dEg^rag/dt in eq. [H] and dE vm /dt in eq. [T5]) and the 
energy dissipation rate (dE tuI b/dt in eq. [24]). In steady 
turbulence, 



dE, 



turb 



at 



= C e ff 



dE, 



g, drag 



dE v 



dt 



dt 



(26) 



where the efficiency factor C c fF represents the fraction of 
the released gravitational energy that is transferred to 
turbulence. Comparison with the numerical simulation 
by JHK06 suggests that C cr f = 0.19 (see fc|5.2l below). and 
we adopt this value in this paper. From equations ([8]), 
([12"]) . (fl8|) . and ([24]) . the viscosity parameter is 



a = 2C cS fiT s 



-'d.turb 



(27) 



Note that the right- hand- side of equation (p?7|) is a func- 
tion of a through / m id or /3 in the "surface densities" (see 
eqs. [13], [IH], [10], [IS]). Before obtaining the exact so- 
lution for a numerically, in the following two subsections, 
approximate solutions are derived. 




mid 



Fig. 1. — "Surface densities" E^dragj ^d,vis> an d S^t^t, in the 
small particle limit T 3 <C 1. The "surface densities" are normalized 
by the dust surface density E^, and are plotted as functions of the 
midplane dust-to- gas ratio, / m id- 



4.1. Turbulence Strength in the Small-Particle Limit 

In equation (l2"7) . the "surface densities", S^drag, 
Eci.vis, and S^.turb represent the effective density of the 
dust that contributes to the liberation of the gravita- 
tional energy or dissipation in the turbulence. The "sur- 
face densities" are functions of the stopping time T s and 
the midplane dust-to-gas ratio /mid- For small dust par- 
ticles (T s <C 1), the "surface densities" in equations (fl4| . 
(|2"0)) , (or approximated eqs. [13], [IS]), and ([21]) de- 
pend on T s only through / mid ~ Z disk (T 1 s /a) 1/2 . Figure 
[T] shows how the "surface densities" vary with / m id for 
small particles (T s <C 1). When the midplane dust-to- 
gas ratio is much smaller than unity, the liberation of 
the gravitational energy mainly comes from the dust ac- 
creting due to individual drag, while for / m id ^> 1, the 
collective drag dominates the energy liberation. For the 
energy dissipated in turbulence, if / m id <C 1, the "surface 
density" E^turb is the column density of the gas in the 
dust layer and is inversely proportional to / m id- If the 
dust layer is thin enough that f mi d 3> 1, then Sd.turb is 
simply given by the dust surface density E^. 

Consider how the turbulent parameter a given by 
equation (|27j) depends on the stopping time T s and the 
disk metallicity Zdisk- In the following discussion, we as- 
sume small particles, such that T s « 1. As discussed 
above, for T s <C 1, the "surface densities" depend on T s 
only through / mid ~ Z disk (T s /a) 1/2 , and consequently 
equation (|2T[) determines the value of a/T s . Thus, in the 
limit of T s <C 1 , it is seen that a cx T s . We derive an ap- 
proximate expression of a for the small and large limits 
of /mid- For / m id -C 1, Figure Q] (and eqs. [15] . [20] . and 

[25]) Show that Srf,drag = Sd, ^d.vis ~ /midS di drag < 

S d ,drag, and Sd.tuxb = 2S d /(0r/ m id)- Thus, the ra- 
tio of the "surface densities" in equation ([27]) reduces 

tO (Ed,drag + S £ j iV i s )/S £ ; i turb =V / ^'/mid/2. For / m id ^ 1) 

though S^^rag is smaller than £d iV i s , it still makes a con- 
tribution to the energy liberation. We fit the functional 
form of (Sd.drag + S t ;,vi S )/S cii turb for 10 2 < /mid < 10 4 by 
a power law form, (X^drag + £<i,vis)/£d,turb ~ ^d,of~f d , 
where 5 = 0.94, and T*d,o — 0.71. Using the above ex- 
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disk 



Fig. 2. — The turbulent viscosity parameter a for various values 
of the stopping time T s plotted as a function of the disk metallicity 
Zdisk- The solid lines show a derived from our model. The dashed 
lines show a from the S98 model. We adopt an efficiency parameter 
of the energy supply, C e ff = 0.19, for our model, and the critical 
Richardson number Ri = 0.8 for the S98 model. 



pressions, equation ([2"T]) becomes 

2 

(Tr 1/2 C eS fjZ disk ) 3 T s for Z disk < y/C cS fj 



2EoC e ft?7^ disk J T s for Z disk > y/C^rj 

and the midplane dust-to-gas ratio (eq.[6]) is 

z 2 



(28) 



fmid 



for Z d[sk < y/Ccfffj 



The scale height of the dust layer is 



(29) 



hd 



(ir 1/2 C eS fjZ d i sk ) 3 h g for Z disk < ,/Ccsfj 

^SoCeff^difk) 2 6 h g for Zdisk > \/C c sfj 

(30) 

In the above approximate expressions (|2"8)l - (|3U|) . the up- 
per line represents / m id *C 1 and the lower line represents 
/mid 3> 1. The transition of these expressions occurs at 
/mid ~ lj which corresponds to Z d i s k ~ y/Ccsf}, and a 
maximum value of a, 



On 



C c sf]T s 



(31) 



4.2. Turbulence Strength for a Small Dust-to-Gas 
Ratio(f mid < lj 

We derive the approximate expression of the turbulent 
strength a as a function of T s and Z d { sk in the limit of 
small dust-to-gas ratio at the midplane, / m i d <C 1, but 
without assuming T s « 1. For / m ; d <C 1, the "surface 
densities" (eqs. [14], [20], and [25]) reduce to Hd.drag = 

2^d{T s + 2)/{T s + 1) , Sd,vis ~ /midS^drag <C E<J, drag, 

and £d,turb = 2Ed/(v / ^"/mid)- Thus, the approximate 
expression of equation (j2"7)l becomes 



-^-Coff^-^disk 



1 + 2T S + 2 
1 + T S (T s 2 + 1) 2 



(32) 
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Fig. 3. — The midplane dust-to- gas ratio / m id- The solid lines 
are calculated for various values of the stopping time T s using our 
model. The dashed line is taken from the S98 model. 

The midplane dust-to-gas ratio and the scale height of 
the dust layer are, respectively, 



/mid — 



2 ^disk 



1 + 2T S (T 2 + 1 



2-| 3 



V^C cff f, 1 + T s T? + 2 



and 



hd = 



— - o ff?7Z d i s k 



r 



1 + 2T S (T s 2 + l) 2 



(33) 



(34) 



4.3. Turbulence Strength and Dust Layer Thickness 

The properties of the turbulent viscosity parameter a 
for T s <C 1 or / mid -c 1 are described in the last two sub- 
sections. For general T s and / m id, equation ([2"T]) must 
be solved numerically. We use the van Wijngaarden- 
Dekker-Brent method (Press et al. 1986) to obtain a 
from equation (|27|) . In Figure [2] the solid lines repre- 
sent the variation of a with the disk metallicity Z d i sk 
for various values of the particle stopping time T s . It is 
seen that, for small dust particles (T s < 1), a is pro- 
portional to T s , as discussed in £14.11 The maximum 
C c sf]Ts appears at Z disk ~ \/C c sf) ~ 10~ 2 (and 
1) for T s < 1. The turbulence weakens as Z d { sk 



/mid 



deviates from \fC c sf] (a cx Z^{^. for Z disk < y/C c sfj, and 
a cx Z^sj. 8 for Z d - lsk > s/Ccsf}, see eq. [28] ) . For particles 
of T s > 1, a peaks at a larger Z d i sk (or / m i d ). 

Figure [3] shows the variation in the midplane dust-to- 
gas ratio / m ; d with the disk metallicity Z d - lsk . For small 
particles (T s < 1), a is proportional to T s , and the mid- 
plane dust-to-gas ratio / m i d cx (T^/a) 1 / 2 does not de- 
pend on T s . For a small dust-to-gas ratio (/ m id <C 1), 



/mid oc Z^ k , and for / mid > 1, / mid cx 
if the disk metallicity is larger than 0.2, then the mid- 
plane dust-to-gas ratio becomes as large as 100, which is 
high enough to induce a gravitational instability in the 
dust layer. This result is consistent with the previous 
works by, e.g., S98. Figure 0] shows that the dust layer 
thickness hd is of the order of 10~ 3 — 10~ 2 times the 
thickness of the gas disk. 

5. COMPARISONS WITH PREVIOUS STUDIES 

In this section, the dust layer model that is based on 
the energetics of dust accretion is compared with results 



Z]£ y . Note that 
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Fig. 4. — The scale height of the dust layer normalized by 
the gas scale height h g . The solid lines are calculated for various 
values of the stopping time T s by our model. The dashed li ne i s 
taken from the S98 model, which is calculated using equation l|38jl . 



neglected. The dust density distribution is given by 

i 



Pa 



Rirjh* 



1 



1 + /mid 



1 



(36) 



The half-thickness of the dust layer, z<j, at which the dust 
density becomes zero, is given by 



Z d = \ Rifj 1 - (1 + /mid)" 



(37) 



To compare with the scale height of Gaussian distribu- 
tion, hd, we define the scale height of the dust distribu- 
tion as the vertical dispersion of the dust particles, 



h 2 



Q d z 2 p d dz 
Q d p d dz 



(38) 



The turbulent parameter asek is estimated from equation 
Q, using T s < 1, as 



from previous work. Previous studies have analyzed the 
detailed physics of the dust layer, including the onset 
of KH or streaming instabilities, and several stabilizing 
effects, such as Keplerian shear, using both analytical 
and numerical approaches. It is of interest to determine 
which properties of the dust layer are reproduced by our 
model and which are missing. 

5.1. Comparison with Previous Analytical Studies 
5.1.1. Comparison with Sekiya (1998) 

S98 analytically solved the density structure of the dust 
layer, assuming that the layer structure was adjusted to 
be marginally unstable to KH instability. To obtain the 
density structure, S98 assumed that the Richardson num- 
ber remained close to the critical value for instability 
throughout the dust layer. S98 considered small dust 
particles that were tightly coupled to the gas (T s < 1), 
and treated the gas and the dust as a single fluid. Then, 
the density structure of the dust layer was calculated for 
various values of the disk metallicity Z d i sk . In S98, the 
dust layer structure was determined by the argument for 
the stability of a stratified fluid, not by the balance be- 
tween sedimentation and diffusion of the dust. However, 
it is possible to interpret the result of S98 as follows: 
The turbulent strength in the dust layer is adjusted such 
that the turbulent diffusion of the dust maintains the 
marginally unstable density structure. Using this inter- 
pretation, we calculate the effective value of the turbulent 
diffusion parameter a as a function of Zdisk f° r T s <ti 1. 

From equation (22) of S98, the midplane dust-to-gas 
ratio, /mid, is related to the disk metallicity, Zdisk, as 



Z, 



disk " 




1 



1 + /mid 



(1 + /mid) 



1 + / 



mid 



aSok 



:,Sck 



(39) 



For / m id > 1, equation ([3"T]) reduces to z d « (Ri^) 1 / 2 /^. 
Substituting this expression into hd,Sek ~ z d of equation 
([391) gives 

a S ck « RWs for /mid > 1 , (40) 

which can be compared with equation (|31|) . For /mid 
1, equation (|35|) is approximated in the lowest order of 
/mid as 



■^disk — 



fn 



4 Rifj 

Thus, the midplane dust-to-gas ratio is 

9 7T 

16Rh7 



7 2 

"^disk 



(41) 



(42) 



For /mid "C 1 and z -C h g , the dust density distribution 
(eq.[36]) is 



pd 
Pg 



fn 



2Ri77/, 



mid 



(43) 



To the order of (z/hd) 2 , this distribution is approximated 
by the Gaussian distribution / m id exp[— z 2 /(2h 2 d gck )] ps 
/mid[l - ^ 2 /(2/i| S ek)]> and its scale height h d ^ Sck is 



1/2 ry 

— 7T 1 Z disk 7] 



hd,Sek ~ y/Rivfmidhg 

The turbulent parameter asek of equation 
"3Ri j 



(44) 



is 



asck 



(35) 



— ir 1/2 Z disk rj 



T s for /mid « 1 • (45) 



where the Richardson number Ri is constant throughout 
the dust layer, and the self-gravity of the dust layer is 



Comparing this expression to equation (|28[) provides a 
relationship between the energy supply efficiency to tur- 
bulence, C e ff, and the critical Richardson number, Ri, as 
C e ff = (3/4)Ri. In the above discussion on asek, we used 
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size (cm) 

Fig. 5. — Comparison with the results in Michikoshi & Inutsuka 
(2006). The solid line shows the relationship between the midplane 
dust-to-gas ratio / m i<j and the particle radius a. The dashed line 
shows the boundary at which the growth rate of the KH instability 
is io = 0.1 (above the dashed line, u) > 0.1). We use the same 
disk model that is adopted in Michikoshi & Inutsuka (2006). The 
dashed line is read from Fig. 16 in Michikoshi & Inutsuka (2006). 

the dust layer thickness /id.Sck derived from the Gaus- 
sian fit. However, the thickness defined by the vertical 
dispersion (eq.[5H]) in the small / m i<j limit is 0.63 times 
thinner than that defined by the Gaussian fit. Thus, 
the energy supply efficiency is 0.4 — 0.63 2 times smaller 
than the above estimate and thus C c g « 0.3RL The re- 
lationship between C e ff and Ri can also be derived by 
comparing the expressions for / m id (eqs. [29] and [42] ) 
as C c ff = 16/(97r 3 / 2 )Ri ~ 0.3RL We compare our re- 
sult with the numerical simulation developed by JHK06, 
which will be discussed in £15.21 The comparison with 
this simulation suggests a slightly smaller energy supply 
efficiency, 

Ccff ~ 0.24Ri . (46) 

The turbulent parameter asek given by equation (|39p is 
numerically calculated and is plotted using dashed lines 
in Figure [5] The critical Richardson number Ri = 0.8 is 
adopted to fit the simulation results provided by JHK06. 
For Zdisk ^ 1 an d T s < 1, the solid and dashed lines 
agree with each other very well. This is expected from 
the above discussion that the dependence of a and asok 
on Zdisk and T s for small / m id are the same (eqs. [35] and 
[4"5"]). For Zdisk ^ 10~ 2 , a deviation between the solid 
and dashed lines appears, and the difference increases 
with i?disk- We note that for Zdisk 0.1, the S98 density 
distribution shows a cusp at the midplane, which is a 
significant difference from the Gaussian distribution of 
our model. In Figure [3J the midplane dust-to-gas ratio 
of the S98 model is plotted with a dashed line. For small 
Zdisk < 0.05, the result from S98 is consistent with our 
result (the dashed line coincides with the solid lines for 
T s < 1). For Zdisk ~ 0.05, the midplane dust-to-gas 
ratio from the S98 model is larger than our results. Both 
results are still qualitatively consistent, indicating that, 
for Zdisk <^ 0.1, the midplane dust-to-gas ratio becomes 
large enough for gravitational instability (/ m id ~ 100). 

5.1.2. Comparison with Michikoshi & Inutsuka (2006) 

Michikoshi & Inutsuka (2006) analyzed the growth rate 
of the KH instability of a dust layer, taking into account 
the relative motion and the friction between the dust and 




disk 



Fig. 6. — The turbulent viscosity parameter a estimated from 
our model, compared with the simulation described in JHK06. The 
solid line is calculated from our model and the dashed line is taken 
from the S98 model. The squares are the simulation results, St in 
Table 2 of JHK06. The stopping time T s = 0.1 is adopted both for 
our model and the simulation. 

the gas. Their formulation does not assume tight dust- 
gas coupling, and thus it can be applied to scenarios in- 
volving large dust particles (T s ^ 1), while the vertical 
gravity, Coriolis force, and Keplerian shear are neglected. 
The initial velocity gradient in the vertical direction is 
caused by the dust inertia, using the formula provided in 
NSH86. The growth rate of instability has been derived 
for a wide range of dust sizes a and midplane dust-to-gas 
ratios / m id, and is summarized in their Figure 16. They 
argued that if the effect of Keplerian shear is taken into 
account, the line corresponding to the growth rate u) = 1 
in the a-/ m ;d plane would be the boundary between the 
stable and unstable configurations. This marginally un- 
stable structure of the dust layer presented by Michikoshi 
& Inutsuka (2006) can be considered an extension of the 
S98 model for general sizes of the dust particles, and 
provides a reference to be compared with our model. 

In Figure [5] we compare our model in the a-/ m id plane 
with the line at which the growth rate u> has a con- 
stant critical value. We adopt the critical growth rate as 
to = 0.1, which is smaller than the value Michikoshi & In- 
utsuka (2006) suggested. The qualitative behavior does 
not differ between uj = 0.1 and lj = 1. Our model well 
reproduces the result of Michikoshi & Inutsuka (2006). 
The agreement of our model with the results of S98 and 
of Michikoshi & Inutsuka (2006) suggests that the onset 
of KH instability is controlled by the energy supply due 
to dust accretion over a wide range of particle sizes, a, or 
stopping times, T s . It must be noted, however, that the 
initial state assumed in Michikoshi & Inutsuka (2006) is 
such that the vertical velocity shear appears only in the 
dust layer, neglecting the velocity shear in the Ekman- 
like boundary layer that may appear between the dust 
layer and the upper gas layer. This is the same assump- 
tion that we adopt. As discussed in ^3.21 this assumption 
is appropriate for small dust particles (T s < 1), while for 
large particles (T s ^> 1) the vertical velocity profile of 
the thick turbulent boundary layer may quickly deviate 
from that given by NSH86. Though the agreement with 
the result by Michikoshi & Inutsuka (2006) shows the 
applicability of our model to KH instability in the dust 
layer for general values of T s , we should keep in mind the 
limitation of the models mentioned above. 



10 



Takeuchi et al. 




0.02 0.04 0.06 0.08 0.1 

Z disk 



Fig. 7. — Midplane dust-to-gas ratio / m id derived from analytical 
models and from the simulation described in JHK06. The solid line 
is calculated from our model and the dashed line is taken from the 
S98 model. The squares are the simulation results, which are read 
from Figure 13 of JHK06. The stopping time is T s = 0.1. 

5.2. Comparison with Previous Numerical Simulations 
5.2.1. Comparison with Johansen et al. (2006) 

Numerical simulations of KH instability induced by 
formation of a dust layer were presented in JHK06. Their 
two-dimensional simulation on the 0z-plane solved for 
both the gas and dust motions. They obtained a quasi- 
steady or oscillating density distribution of the dust in 
which dust settling and turbulent diffusion balanced each 
other. The diffusion coefficient in the vertical direc- 
tion, S t , was measured from the scale height of the dust 
layer and was summarized in their Table 2. In Fig- 
ure [6l the turbulent viscosity parameter a of equation 
(|27|) (the solid line) is compared with the diffusion coef- 
ficient measured from the simulation (the squares). To 
fit the simulation result, the efficiency parameter of the 
energy supply to turbulence, C c g, is set to 0.19. We 
also fit the turbulent diffusion parameter for the density 
distribution of S98, ctsck, to the simulation. We adopt 
the critical Richardson number as Ri = 0.8 for fitting. 
This is slightly smaller than the value JHK06 suggested 
(Ri = 1.0). It is seen that the numerical simulation is well 
explained both by our model and the S98 model. This 
means that the dust layer is maintained so as to provide a 
constant Richardson number, and in the parameter range 
that JHK06 surveyed (0.01 < Z disk < 0.1), this critical 
Richardson number does not vary significantly with the 
disk metallicity .Zdisk- (The recent result by Lee et al. 
(2010), in which they argue that Ri is proportional to 
•Zdisk if the Keplerian shear is taken into account, will be 
discussed later in this subsection.) Because the efficiency 
parameter C c g in our model is proportional to the crit- 
ical Richardson number, as shown in equation (|46p. the 
efficiency parameter C e g is also expected to be constant 
for 0.0K Z disk < 0.1. 

A comparison of the midplane dust-to-gas ratio be- 
tween our model and the simulation is shown in Figure 
[7J The simulation is well fitted by our model (the solid 
line) and also by the S98 model (the dashed line). For 
-^disk = 0.1, the S98 model predicts larger / m ;d than the 
simulation. This may be because the simulation does not 
have enough resolution to resolve the density structure 
around the midplane at which the S98 model expects a 
rapid increase in the dust density. If this is the case, 
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Fig. 8. — The scale heights of the dust layer estimated from 
our model and from the JHK06 simulation are plotted versus the 
stopping time T s . The solid line is calculated assuming that the 
turnover time of the largest eddies is the Keplerian time (T e = 1). 
The dashed line includes the variable T e calculated by equation 
l|48p. The squares are the simulation results, (z 2 ) 1 / 2 /H in Table 2 
of JHK06. The disk metallicity is Z disk = 10 -2 . 

it is difficult to judge whether our model or the model 
presented by S98 is the best fit with the simulation. 

Figure [H] shows the variation in dust layer thickness 
with T s from our model and compares this variation with 
the simulation results. In plotting Figure |6] in the last 
paragraph, the turbulent viscosity parameter a was cal- 
culated from equation Q using the dust layer thickness 
hd measured from the simulation. However, it is not 
clear if equation (@| is still valid for large T s because this 
equation assumes that the turnover time of the largest 
eddies is equal to the Keplerian time (YL07), and as dis- 
cussed below this assumption may not be appropriate 
for large T s . Thus, in Figure El the dust layer thick- 
ness is plotted directly without transferring to a. In 
our model, the dust layer thickness hd is constant for 
T s < 1 and decreases slowly with T s for T s > 1 (the 
solid line), although in the simulation it decreases more 
rapidly with T s . The discrepancy between our model 
and the simulation is apparent for T s = 1. The differ- 
ence is as large as a factor of 3, but it causes an order 
of magnitude discrepancy in a a /i^. This discrepancy 
suggests that our model predicts turbulent diffusion that 
is too large for T s > 1. In fact, for large T s in our model, 
the dust layer thickness becomes smaller than the size 
of the largest eddies, contradicting our assumption that 
the turbulent layer coincides with the dust layer (see, 
however, the discussion in H3.2l on the validity of this as- 
sumption, and see also eq. [56] of YL07 for a possible 
physical reason for ^ ff , c ddy > hd)- From equation ((4]) and 
^s.oddy = OL l l 2 h g , the condition for the dust layer thick- 
ness to be larger than the largest eddy size (hd > l g , e ddy) 
is T s < 1/V2. Thus, for T s > l/y/2, our model results 
in a dust layer that is too thick (or an eddy size that 
is too small). One possible remedy for this situation is 
to remove the assumption that the turnover time of the 
largest eddies is equal to the Keplerian time. By intro- 
ducing a parameter £ (0 < £ < 1/2), the largest eddy size 
and the velocity are expressed as / s , e ddy = cr 1 ' 2 ) + ^/i 9 and 
w ff ,eddy = or- l ^~^c a . The non-dimensional turnover time 
of the largest eddies is T e = (? 9 ,cddy/wg,eddy)^K — 
The dust layer thickness for T e ^ 1 is provided by equa- 
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tion (21) of YL07, 

^vfC+S^' <47) 

For T s > 1/2, we impose the condition that the largest 

eddy size l g , e dd y = a^^ 2 ' + ^h g is equal to the above hd- 
This condition determines the turnover time by the equa- 
tion 

T s T e (T s +T e + T s T e 2 )-T s -T e =0 . (48) 

Using this T e (or £) , the turbulent viscosity parameter a 
and the dust layer thickness hd are recalculated and plot- 
ted with the dashed line in Figure [SJ As expected, the 
introduction of a new parameter T e suppresses the largest 
eddy size and the turbulent diffusion of dust particles for 
T s > l/v2j improving the comparison with the simula- 
tion result. However, we need to determine whether the 
turnover time of the largest eddies, Tefl^ 1 — o^f^ 1 , in 
the simulation presented in JHK06 for T s = 1 is actu- 
ally smaller than the Keplerian time, as equation (|48|) 
predicts. Cuzzi et al. (1993) argue that, based on labo- 
ratory measurements, the eddy turnover time in an Ek- 
man layer would be smaller than the Keplerian time by 
a factor 20 — 80. The turbulent layer possibly behaves 
as an Ekman layer for T 8 > 1, as discussed in i)3.21 We 
are currently performing numerical simulations using the 
same conditions as JHK06 to investigate the turbulence 
for T s > 1 in more detail (Ishitsu et al., in preparation). 

Lee et al. (2010) performed a three-dimensional nu- 
merical simulation of the onset of KH instability. They 
solved simplified equations in which the dust and the gas 
were treated as a single fluid, but they included the ef- 
fect of the Keplerian shear in the radial direction. They 
found that the radial shear stabilizes the KH instability, 
and the critical Richardson number for instability is not 
always the standard value, 0.25, but can be much smaller 
if the stabilizing effect of the radial shear is significant 
(sec also Ishitsu & Sekiya 2003). Equation (32) in Lee et 
al. (2010) shows that the ratio of the stabilizing effect by 
the radial shear to the destabilizing effect by the verti- 
cal shear is proportional to Ri(l + / m id)//mid ~ Ri//mid 
for / m id *C 1, and thus the critical Richardson number 
should scale as Ri oc f m id- In our model, we assume 
that the energy supply efficiency C c ff, which is propor- 
tional to the critical Richardson number (eq.|46j). is con- 
stant. However, the simulation in Lee et al. (2010) 
suggests that C c g should also be proportional to / m id 
(even for / m id > 1). If this is the case, dependence of 
/ m id on Zdi s k would be milder than that shown in Fig- 
ure [3] (/mid oc Z\r k for /mid < 1 and / mid oc Z^ for 
/mid ^1). The effect of the radial shear was not included 
in the simulation presented in JHK06, with which we 
compared our model in detail, because their simulation 
was two-dimensional in the #z-plane. Extending JHK06 
to three dimensions and including the radial shear ef- 
fect are crucial to determining how the energy supply 
efficiency C e s behaves as / m id varies. 

5.2.2. Comparison with Bai & Stone (2010) 

Bai & Stone (2010) performed a three-dimensional sim- 
ulation, focusing on investigating turbulence induced by 
streaming instability. They found that the streaming in- 
stability induced turbulence before the KH instability set 
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Fig. 9. — The turbulent viscosity parameter a estimated from 
our model for T s = 10 _1 and 1 is compared with the simulation 
presented in Bai & Stone (2010). The solid line is calculated as- 
suming the turnover time of the largest eddies is the Keplerian 
time (T e = 1). The dashed line includes the variable T e calculated 
from equation 1148 II . (a) Comparison with the vertical diffusion 
coefficient, D 9iZ (3D), measured in the simulation (Table 2 of Bai 
&; Stone 2010). The squares and circles are the simulation re- 
sults, connected by lines that indicate which runs have the same 
size distribution of dust particles. The label of run Rxj/ means 
the particle size distribution is such that the stopping time ranges 
T s e [10- x ,W- y ], e.g., in run R30, T s e [10^,10°]. In each 
run, the disk metallicity Z is 0.01 and 0.03. To compare the re- 
sults with our model, we assume that only the largest particles 
contribute to turbulence, and that the effective metallicity is esti- 
mated by Z d ; sk = Z/A?typ e , where N typc = 3 (R10 and R21) or 7 
(R30 and R41) is the number of particle species. We compare our 
model of T a = 1 with runs R10 and R30 (plotted in blue), and the 
model of T s = 10 _1 with runs R41 and R21 (plotted in red), (b) 
Comparison with the radial diffusion coefficient, D x , measured in 
the 3D runs of the simulation. We assume that the diffusion co- 
efficient of the gas is represented by that of the smallest particles 
in the simulation. The values of D x for the smallest particles are 
read from Fig.9 of Bai & Stone (2010). 

in. They measured the turbulent diffusion coefficient. 
However, it is difficult to compare our model directly 
with the results in Bai & Stone (2010), because their 
simulation includes particles of several sizes (3-7 species) 
while our model considers only single-sized particles. In 
the simulation, it was reported that only large particles 
were responsible for inducing turbulence. Our model also 
shows that the energy liberation per unit mass of the 
dust is higher for larger particles (it is proportional to 
T s for T s < 1). In order to compare results, we assume 
that in the simulation, the turbulence is induced only 
by the largest particles (i.e., particles of largest T s ). For 
example, in the R41 run (in which the stopping time of 
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the particles ranges from 10 -4 to 10 _1 ), we assume that 
only T s = 0.1 particles are responsible for turbulence. 
We then compare the simulation result with our model 
of T s = 0.1. In the simulation, each species has the same 
amount of mass. Because we consider the largest parti- 
cles only, the total amount of the dust participating in 
driving turbulence is Ed/-/Vtyp C , where N type is the num- 
ber of particle species in the simulation. For example, 
in the R41 run, Nt ypc = 7, and we compare the simula- 
tion with a disk metallicity Z = 0.01 with our model of 
Z disk = Z/N typc = 1.43 x 10~ 3 . 

Figure Eh shows a comparison of the diffusion coeffi- 
cient obtained from the simulation (D g z (3D) in Table 2 
of Bai & Stone (2010)) with our model (calculated with 
the parameter f\ = 0.05 2 , which was adopted by Bai & 
Stone (2010)). Although our model of T s = 0.1 agrees 
with simulations R41 and R21 (plotted in red), we note 
a qualitative discrepancy between the T s = 1 model and 
simulations R30 and R10 (plotted in blue). The sim- 
ulations indicate that the diffusion coefficient decreases 
with the disk metallicity, and that its value for Z = 0.03 
is about an order of magnitude smaller than the value 
from our model. This discrepancy cannot be resolved, 
even by varying the turnover time of the largest eddies 
T e (the dashed line). One possible cause for the inconsis- 
tency is particle clumping and concentration in turbulent 
eddies, which are not included in our model. The simu- 
lation shows strong particle clumping when Z — 0.03 in 
the R10 run and also temporal clumping for Z = 0.03 
in the R30 run. Such clumping of particles in turbulent 
eddies may suppress diffusion of particles compared to 
the no-clumping cases of Z = 0.01 and could be a cause 
of a decrease in a when Z is increased in the simulation. 
Even if the largest particles (of T s ~ 1) concentrate in 
clumps, the smallest particles (of T s < 0.1) do not clump, 
and continue to follow the turbulent diffusion of the gas 
(Fig. 7 of Bai & Stone 2010). In Figure©, the tur- 
bulent diffusion coefficient in the "radial direction" of 
the smallest particles in the simulation, D x , is compared 
with the "vertical" diffusion coefficient a in our model. 
Note that we compare diffusion coefficients in the differ- 
ent directions. Since the smallest particles spread out to 
high altitudes where turbulence is weak, it is difficult to 
measure the vertical diffusion coefficient for the small- 
est particles in the simulation. In Figure M>, though we 
still see a discrepancy compared with the R30 run, our 
model appears more consistent with the simulation re- 
sults, suggesting that our model properly predicts the 
"gas" diffusion coefficient. 

6. DISCUSSION 

6.1. The Radial Drift Velocity and Collision Velocity of 
Dust Particles 

If the dust-to-gas ratio in the dust layer were larger 
than unity, the radial drift velocity of dust particles 
would be lower than the value that the particles would 
have in a gas-rich environment because the gas drag force 
could not accelerate sufficiently against the large iner- 
tia of the dust. This e ffect w as pointed out by NSH86 
and is seen in equation (IA14[) for the radial velocity due 
to individual drag. If the dust-to-gas ratio were much 
larger than unity, dust accretion would be caused by col- 
lective drag exerted from the slower-orbiting upper gas 
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Fig. 10. — The radial drift velocity, Vd r (solid line), and the 
relative velocity of particles due to turbulence, Av±2 (dashed line). 
The contributions of individual drag, i)j n d,r> an d collective drag, 
v coi r , to the radial drift velocity are plotted with dotted lines. 

layer (Weidenschilling 2003). The radial drift velocity of 
the dust is thus a function of the turbulence strength, 
a. For weaker turbulence, the dust-to-gas ratio in the 
dust layer is higher, and the individual drag is weaker. 
The collective drag is also weaker at smaller a because 
the Reynolds stress Pg z is proportional to a. The rela- 
tive velocity (or collision velocity) of the dust particles is 
also a function of a. The radial drift and collision veloci- 
ties are the important factors in the dust growth process. 
In the previous sections, the turbulence strength a and 
the dust layer thickness hd have been determined self- 
consistcntly. Using these results, the radial drift velocity 
and collision velocity of the dust particles are estimated. 

The radial drift velocity is calculated separately for the 
components due to individual drag and due to collective 
drag. For each component, the radial drift velocity is 
averaged in the vertical direction. First, the averaged 
value of the radial drift vel ocity due to individual drag 
is calculated from equation (|A14|) . 

1 



"ind.r = — - / PdVd,rdz 



-2i]v K T s — 

Ed J_ c 



1 



Pd 



T s 2 + /3 : 



-dz 



(49) 



The radial drift velocity due to collective drag is cal- 
culated from the vertically-averaged angular momentum 
loss of the dust component. Integrating equation (|C4|) 
gives, 



dL d 
dt 



did 

-no 9t 



dz = -2j7«|T s £ d) , 



(50) 



where Ed, vis is given by equation (|C10|) . This angular 
momentum loss causes a radial drift velocity v co \ r given 

by 

2 dLd . m Erf vis /-..n. 
Vco\,r = =--7T = ~^V v kT s — . (51) 

v K 2^d dt E d 

The total radial drift velocity is Vd, r = Vi n d,r + v co i.r, and 
is shown in Figure [TU] as a solid line, for the case in which 
the particle size is chosen to maximize the radial velocity 
(T s = 1). In plotting this figure, we adopt the model 
parameters at 1 AU of the minimum-mass-solar-nebula 
model of Hayashi (1981): h g /r = 0.0326, rj = 1.80xl0~ 3 , 
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Fig. 11. — The midplane dust-to-gas ratio / m id is plotted against 
fj for various values of the disk metallicity Zdisk- The top axis 
shows the corresponding radii in the gas disk model by Hayashi 
(1981), in which the disk temperature is T = 278(r/AU)- 1 / 2 K. 
The dust particles are assumed to be small (T a -C 1) such that 
/mid is independent of T s . 

and fj = 2.92 x 1CP 3 . For a small disk metallicity Zdisk) 
the radial drift velocity is as large as 50 m s" 1 , and it 
decreases with Zdisk- For Zdisk > 0.08, the radial velocity 
is dominated by collective drag (see the dotted lines) as 
pointed out by Weidenschilling (2003). For such large 
Zdisk, the radial drift velocity due to collective drag also 
decreases with Z^^, and then it becomes as small as 1 m 
s _1 for Zdisk = 0.2. The radial drift velocity is strongly 
suppressed for a sufficiently massive dust layer. 

The relative velocity of dust particles due to turbu- 
lence, A i/i2, is calculated for T s = 1 by substituting 
Sti = 1 and St2 = into equation (29) in Ormel & 
Cuzzi (2007), and is shown in Figure [10] as a dashed line. 
The collision velocity is estimated by the larger of Av\2 
and Vd, r - For small disk mctallicities Zdiskj the collision 
velocity is dominated by the radial drift and is as large as 
50 m s _1 , while for large Zdisk > 0.02, it is dominated by 
turbulence. The maximum value of the collision velocity 
due to turbulence is about 30 m s _1 , and it decreases 
with Zdisk for Zdisk > 0.03. Thus, if the dust parti- 
cles could survive collisions of 30 m s _1 , as suggested by 
the numerical simulation of collisions of dust aggregates 
(Wada et al. 2010) and Zdisk }t 0.03, the dust particles 
would be able to grow without being reduced to small 
fragments. 

6.2. Radial Dependence of the Midplane Dust-to- Gas 

Ratio 

As shown in jgj the turbulence strength and the dust 
layer structure depend on properties of the gas disk only 
through fj. For larger fj, the accretion velocity of the dust 
is faster, and consequently the dust layer is thicker due to 
stronger turbulence. Figure [Til shows how the midplane 
dust-to-gas ratio varies with fj. The figure is plotted in 
the limit of T s <C 1. For such small particles, / m id is in- 
dependent of T s (see Fig. [3]) . The midplane dust-to-gas 
ratio / m id decreases with fj as shown in equation (|29p . For 
a gas disk model with a power-law temperature profile, 
T oc r~ q , fj behaves as fj oc r 1 ^ 9 . The top axis of Figure 
ITTI indicates corresponding locations in the disk model by 
Hayashi (1981), i.e., T = 278(r/AlfT 1 / 2 K. In a disk with 
the standard value of the disk metallicity Zdisk = 0.01, 
/mid is less than unity in the most part of the disk except 



r < 0.1AU. In disks with Zdisk = 0.3, / m ;d exceeds 100 for 
r < 1AU. We discuss the condition for planetesimal for- 
mation through gravitational instability of the dust layer. 
In a gas disk with a surface density profile E 9 oc r~ p , the 
midplane gas density scales as p g oc r~ p+q / 2 ~ 3 / 2 . The 
condition that the dust density exceeds the Roche density 

(p R OC r~ 3 ) is /mid > p R /p g OC r P-9/ 2 -3/2 ^ r P-7/4 ft 

disk model with q = 1/2. From the lower line of equation 
([29]). / mid decreases as / mid oc r -(i-<?)/(2-<5) x r -0A7 _ If 
p > 1.3, the inner part of the disk is more suitable for 
gravitational instability, and vice versa. 

6.3. Does the Liberated Gravitational Energy Go into 
Turbulence ? 

In this paper, we calculate the liberated gravitational 
energy from dust accretion, assuming some fraction of 
the liberated energy is transferred to turbulence. The 
estimate of the dust accretion rate is based on the for- 
mula for the particle terminal velocity derived in NSH86. 
"The terminal velocity" means that all the liberated en- 
ergy is consumed by gas drag, converting directly into the 
thermal energy of the dust particles and of the surround- 
ing gas molecules. Thus, one may expect that only a 
small fraction (or nothing) of the liberated energy would 
be used for maintaining turbulence. However, a com- 
parison with the simulation of KH instability by JHK06 
shows that the efficiency factor C e g « 0.2 is not negli- 
gibly small. In the following subsection, we discuss the 
validity of using the particle terminal velocity for calcu- 
lating the energy supply rate to turbulence. The energy 
liberation rate from the accreting dust calculated in <J3]is 
compared with the deposit rate of the free energy that is 
the source of several instabilities, such as KH instability 
and streaming instability. 

6.3.1. Kelvin- Helmholtz Instability 

The free energy that induces the KH instability origi- 
nates from the velocity difference between the midplane 
dust layer and the upper gas layer, and is stored as the 
dust particles settle to the midplane. We estimate the 
deposit rate of the free energy during dust sedimenta- 
tion, and show that it has the same order of magnitude 
as the energy liberation rate from the dust accretion to- 
wards the star. Consider two states of dust distribu- 
tion: the initial state, in which the dust particles are 
distributed uniformly in the gas disk, and the final state, 
in which all the dust has settled at the midplane. In 
the initial state, there is no vertical shear in the disk, 
and in the final state, the velocity difference Ave — ?1 v k 
appears between the midplane dust layer and the up- 
per gas layer. The free energy for KH instability is 
AE KH ~ \^ d Av 2 ~ for S d < S 9 . The 

settling timescale is r sc d = (T s 2 + 1)/(T s CI-k), and then 
the energy deposit rate is 

AEkh 1 2 2 r> Ts V (M\ 

~ 2 V ^lf+i^ ' (52) 

which is the same order as the energy liberation rate of 
the accreting dust (eqs. [T2] and [IS])- Hence, the depo- 
sition rate of the free energy for KH instability can be 
estimated by the energy liberation rate of the accreting 
dust. 
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6.3.2. Streaming Instability 

The free energy for streaming instability originates 
from the velocity difference between the dust particles 
and the surrounding gas. When streaming instability 
begins, the velocity difference decreases as the free en- 
ergy is consumed by inducing turbulence. In fact, this 
decrease in the velocity difference can be seen even in 
the linear growth regime. Youdin & Goodman (2005) 
showed in their Figure 6 that the velocity difference be- 
tween the dust and gas decreases (increases) at the lo- 
cations where the particle density increases (decreases). 
The spatially averaged value of the free energy decreases 
as the perturbation grows. Thus, without an energy sup- 
ply, streaming instability would cease. Given a state in 
which the velocity difference between the dust and the 
gas has reduced, the dust particles are no longer in equi- 
librium: the gravity, the centrifugal force, and the gas 
drag force are not in balance. The dust particles are ac- 
celerated and the velocity difference from the gas rises 
again. To estimate the effect of the energy deposition on 
the velocity difference, we consider a state in which the 
dust density is similar to the gas density, pd ~ p g - In 
such a state, streaming instability occurs efficiently with 
a growth time of the order of the Keplerian time (for the 
short wave branch, Youdin & Goodman 2005; Youdin & 
Johansen 2007) . For dust particles of T s < 1 , the termi- 
nal velocities of the dust and of the gas are of the order 
of T s r)v K (eq. [S], [TU]), and thus, the free energy per 
unit area is AE stl - ~ T.JTgrfv 2 ^. This deposition of free 
energy occurs during the acceleration phase of the dust, 
and thus in the stopping time Tsfi^ 1 , and then the energy 

is transferred to turbulence in the growth timescale Q^ 1 
of streaming instability. Thus, the timescale for the en- 
ergy deposit in turbulence is the sum of these timescales, 
and for T s < 1, r str = T^ 1 + fly- 1 ~ Q,^ 1 . The energy 
deposition rate is 

^ ~ V 2 v 2 K n K T*Z d , (53) 

Tstr 

which is smaller by a factor T s than the estimate from 
the dust accretion rate (eq.[l2]). Hence, our estimate, 
derived from the dust accretion rate, is appropriate for 
particles of T s ~ 1. For smaller particles (T s <C 1), the 
energy deposition rate is higher for KH instability than 
for streaming instability, and KH instability is expected 
to operate first. The energy deposition rate for instability 
(cither for KH or streaming instabilities) is estimated 
from the energy liberation rate due to dust accretion. 

7. SUMMARY 

In this paper, we discuss turbulence induced in the 
dust layer. The turbulence strength or the parameter 
a is determined using the energetics of dust accretion 
towards the central star. The key concept is that the 
dust particles reside in a deeper potential than the gas. 
The effective potential, including the gas pressure, is 
—GM(1 — 2rj)/r for the gas, and —GM/r for the dust. 



When angular momentum is transferred from the dust 
to the gas through gas drag, the dust particles lose more 
energy than the gas gains. The excess energy can be used 
for exciting turbulence. If the dust accretion due to gas 
drag is a primary source of energy liberation, i.e., if the 
gas accretion rate due to turbulence is much smaller than 
the dust accretion rate, then the turbulence strength is 
determined by the energy supply rate from the dust ac- 
cretion. This is not the case if the gas disk itself is tur- 
bulent via, e.g., MRI. If the dust layer is composed of 
large particles with stopping time T s 3> 1, then the gas 
accretion may dominate the dust accretion, as discussed 

in m 

We estimate the dust accretion rate using the terminal 
velocity profiles of the dust particles in a laminar disk 
derived by NSH86. The expected turbulence strength 
and corresponding structure of the dust layer from our 
analysis agree with the previous analytical result on the 
marginally KH- unstable dust layer by S98. As our anal- 
ysis does not assume tight coupling of the dust to the 
gas, nor specify the mechanism of instability that induces 
turbulence, it is considered an extension of the analysis 
of S98 to a more general physical situation of the dust 
layer. The results of this paper agree with the results in 
Michikoshi & Inutsuka (2006), which analyzes KH insta- 
bility of the dust layer composed of particles with large 
stopping times (T s > 1), as shown in Figure [5] 

Our analysis shows that, for particles of T s < 1, the 
turbulence strength is smaller than a max ~ C c fifjT Sl 
where C c s = 0.19 is the efficiency of the energy sup- 
ply to turbulence (see Fig. [5] and eq. [SI]). This 
strength reaches a maximum when the disk metallicity is 
-Zdisk ~ y/Cctffj ~ 10~ 2 . Modifying the disk metallicity 
from the standard value, 10 -2 , by any process, results 
in weaker turbulence and a thinner dust layer, and con- 
sequently may accelerate the growth process of the dust 
particles, as pointed out in S98. 

Comparison of our results with previous numerical sim- 
ulations of KH and streaming instabilities by JHK06 and 
Bai & Stone (2010) shows quantitative agreement with 
our analysis for dust particles of T s < 0.1, although there 
may be a qualitative disagreement for T s > 1 particles 
(see Figs. [B]— [5]). Hence, we conclude that turbulence in 
the dust layer is controlled by the energy supply from the 
dust accretion due to gas drag, provided that the dust 
particles are not so large that T s > 1. In such a layer, 
turbulence strength is estimated by the dust accretion 
rate (eq.gTj). 
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APPENDIX 

DUST AND GAS VELOCITIES IN STEADY LAMINAR FLOW 

We present calculations of the dust and gas velocities in steady laminar flow in this appendix. We follow NSH86, but 
extend their calculation to the second order of r\. The equations of motion of the gas and of the dust are, respectively, 

(«-.). <A2) 

where v g and v d are the velocities of the gas and the dust. The radial and azimuthal components of the velocity in 
the cylindrical coordinates (r, 6) are normalized by the Keplerian velocity vk, such as v g . r = v g<r VK, v g< g = v 9i 6Vk- 
We assume that the velocities, v g . r etc., vary with r in the same way as the Keplerian velocity vk, i.e., that the 
non-dimension al velocities , v g , r etc., are constant with r. This assumption is satisfied when pd/pg and 77 are constant 
with r.(see eqs. |Al2] - |A19| below) Then, the radial derivative of the velocity is, for example, 

d 



6r V " = "IF ' (A3) 

and the radial derivat ive of oth er v elocity components has a similar form. In a steady axisymmetric state (d/dt = 
8/86 = 0), equations ([AT]) and (|A"2j) become 

- \&g, r ~ v 2 gfi = -(1 - 27?) - ~(v g ,r " Vd,r) , (A4) 

* Pg J- s 



1- - Pd 1 f . 
2 W = -y g T t " 



Vg,rV g fi = 7^{Vg,e ~ V d fi) , (A5) 



- \v\ r - v\ e = -1 - — (V d , r ~ Vg,r) , (A6) 

l;Vd,rVd,8 = -jr{v d ,o - v g ,e) ■ (A7) 
The non-dimensional velocities are expanded in a power series of 77, 

Vg,r = V g ,r,lV + %r,2?7 2 + 0(l] 3 ) , (A8) 

v g ,e = 1 + v 9t e,iV + v g ,e,2V 2 + 0{r] 3 ) , (A9) 

"d,r = Wd.r,l ? / + Vd.raV 2 + 0{rf) , (A10) 

= 1 + v d ,e,ir] + i d ,e,2V 2 + 0(rf) . (All) 
Substituting the above expressions into equations (|A4j) - (|A7j) yields in the first order of 77, 

Pd 2T S 



,J gfi.\ 



PgTi+P 2 ' 
_T[±£ 



(A12) 
(A13) 



2T 

^ = -jf^ - (AH) 

= - T 2 + p2 » ( Al5 ) 



which are the same as the results of NSH86. In the second order of 77, 



p d Ti(3T 2 + 2f3* + 2!3) 

Vr,2 = q , (A16) 

9 " Pg (T2 + /32) 3 



r 6 6 + 3/?r s 4 + 3/? 2 r s 2 + /? 



4 



(A17) 



^, 2 = - T - (r - +2 ^ , (A18) 



16 



Takeuchi et al. 



Vd.6.2 



(3 (3T S 4 + 2,3 2 T 2 + 3f3T 2 + (3 3 



2(T 2 



n 2 f 



(A19) 



ENERGY LIBERATION RATE DUE TO INDIVIDUAL DRAG 



In this appendix, we describe a more rigorous derivation of the energy liberation rate due to individual drag than 
was provided in £13.11 In a laminar disk, the particle drift velocity Vd, r calcula tion s were presented in NSH86, to the 
first order of rj. Since the liberated energy is of the second order of rj (see eq. |B4| below), we use the particle radial 
velocity, which is calculated to the order of rj 2 in Appendix [3] 



2T„ 



T 2 +/? 2 



T s 3 (T s 2 + 2/ 3) 
(T 2 + /3 2 ) 



2^ 



VK , 



(Bl) 



where the first term corresponds to equation (2.11) in NSH86. The gas drifts in the opposite direction with the velocity 



J 9,r 



Pd 

Pg 



2T S T s 3 (3T 2 + 2/3 2 + 2/?) 2 

V+ 7^73 -V 



V K ■ 



The liberated gravitational energy per unit surface area of the disk is, 

9E dlaK 1 



dt 



(PdgdVd,r+Pg9gVg, r )dz , 



(B2) 



(B3) 



where the factor of 1/2 accounts for the work used for the acceleration (and dec elerati on) of the azimuthal velocity of 
the dust (and of the gas) as their semi-major axes change. Using equation (|Blj) and (|B2|) . the energy liberation rate 
is given by 

= 2ry 2 i; 2 i f!KT ;j £ (; , drag , (B4) 

where the effective "surface density" of the dust is 

S d r exp(-z 2 ) 



T 2 + f3 2 



T 2 



2(T 2 + /3 2 ) 



dz 



(B5) 



ENERGY LIBERATION RATE DUE TO COLLECTIVE DRAG 

In this appendix, we provide a more rigorous derivation of the energy liberation rate due to collective drag than was 
given in £13.21 The orbital velocity of the gas presented in NSH86 is 



V 9fi = 1 ~ 



P + T 2 
P 2 +T 2 ' 



(CI) 



where /3 = (pd + Pg)/ Pg varies with the altitude z. The Oz component of the Reynolds stress Pg z due to the turbulent 
viscosity v of the gas is 

Pez = [Pg + C stT p d )u- ■ 



dz 



(C2) 



We add the factor C s t r to account for the weaker coupling of the dust to the gas for larger dust particles. YL07 has 
shown that the contribution of the dust to the r^-component of the Reynolds stress, P r e, is a factor 1/(T 2 + 1) times 
the gas contribution. From equations (33c) and (B.l) of YL07, it is seen that in both limits of T e <C 1 and T e 1, 

■ d Qj ~ (v' g . r v' g (^jT~ 2 for T s » 1, and (v' dr v' d ^j ~ ( v 'g,r v 'g $) f° r T s <C 1, where the prime denotes velocity 

fluctuations. We assume that a similar relationship holds for the #z-component of the Reynolds stress, Pg z , and thus 
C str is expressed at0 

Cstr ~ 



T 2 + l 



(C3) 



6 Equation (B.l) of YL07 is based on the radial shear effect, and 
its applicability for the 02-component is not very clear. We simply 
assume that P r g and Pg z have similar properties. Note also that 
C s tr includes only the effect of T s , assuming that the dust parti- 
cles act as passive particles in the gas turbulence. This is not the 
case if the local dust-to-gas ratio is larger than unity. Since the 
effect of the inertia of the dust on the turbulence is unclear, we 



simply adopt equati on dC3|l . For / m id > 1, the simple plate drag 
approximation (eq. |19|^ may provide a more accurate esti mate . 
The Ed v i s estimated from the plate drag approximation (eq. 1191 ) 
and from the calculation in this Appendix (eq. ICIOI ) does not 
suggest a big difference at large / m i(j. In the plate drag approxi- 
mation, v ; s oc / — h, while the calculation in this Appendix gives 



c</-.° d M (see Fig. [TJ. 
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This stress transfers angular momentum in the z-direction, and the time derivatives of the angular momentum of the 
dust and of the gas per unit volume and unit time are, respectively, 

dl d C str pd dP Sz 
dt p g + C stl pd dz 

and 

§ = — ~£r r ~7f 1 ( C5 ) 

dt p g + C stI p d oz 

Here, we assume that the viscous torque is distributed to the dust and to the gas with the ratio C stI Pd '■ Pg- The 
corresponding energy change is ded/dt = fl^dld/dt for the dust and de g /dt — Q g dl g /dt for the gasQ. In sum, the 
energy liberation rate per unit area and unit time is 



dE vis f°° ( de d de g 



dt I _ \ dt dt 



OO 

oo 



dz 



= -vkV I Pezi- I ) dz > ( C6 ) 



— OO 



dz \ l + C s tr{pd/Pg) 



where we use P# z (±oo) = 0, and the minus sign is added to give the energy release rate. The integration variable is 
transferred to z = z/(\/2hd), and then using 

9v g ,e - , ~2\ P 2 + 2f3Ts — T 2 

2/ mid 2;exp(-z ) — /rr2 ■ Q2N2 — r]v K , (C7) 



and 



equation (1C6[) becomes 



dz -"—-v > (T 2 +/3 2 ) 

d ( 1 \ 2C str /mid5exp(-z 2 ) 



(C8) 



dz V 1 + Cstiipd/Pg) J [1 + Cstr/mid exp(-Z 2 )] 2 

2r?v 2 K n K T s i: d y la , (C9) 



dEvis 2 2 



where 



dt 

_C stI 1 + 2T S [°° /? 2 + 2/3T 2 - T 2 z~ 2 cxp(-2z 2 ) 

d,ViS ~ ^ TTlT i-oo (T s 2 +/3 2 ) 2 [l + C str / mid exp(-z 2 )] rfZ ■ (C10) 

ENERGY DISSIPATION RATE OF TURBULENT DUST MOTION 
The energy dissipation rate of the gas in turbulence is 



U 3,eddy 



r ff:Cddy 



(Dl) 



where u g cddy and r ff]eddy are, respectively, the velocity and the turnover time of the largest eddies. Similarly, for the 
dust, 

e d = . (D2) 

7"(i,eddy 

The eddy velocity of the dust, w^eddy, is estimated from equation (20) of YL07 as 

2 U l eddy 

<cddy = 1+Ta £l V +TaTe ' ( D3 ) 

where T e = r gi eddy^K is the non-dimensional turnover time of the gas turbulence. The turnover time of the dust 
turbulence, r^eddyj is estimated as the larger value of T 9 cddy and T st0 p, i.e., 

T~d,eddy = max(r 9iCddy ,r stop ) . (D4) 

From equations (|DT|) - (lD4t . 

£d — C cnc £ g , (D5) 

7 We assume that the rotational velocities of the dust and the ,< , . , . , , , , ,. ... n , . , , ■ . 

„ ,„ ,-i i j . ji ,. r ii the dust-dominant layer rotating with S2k) and the gas-dominant 

gas are S!k and S2 respectively, neglecting the modification ol the , , . .. .,, A % ^ ^, t ~ -. ,v' ,. , 

° d .. , mi ■ . . . . ,.„ , . r layer rotating with sZ„). tor tm\d < 1, the energy liberation due 

rotational velocity due to gas drag. This assumption is justified if tQ ^cctivc | is ne 9 g ' lected comp ared to that due to individual 

/mid > 1 and it the angular momentum exchange occurs between drag (see Fig [TJ 
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where C ene is approximately 



T S T { + 1 
T 3 2 (Tr 2 + l) 



for T s 
for T„ 



(D6) 
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